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We report on the application of chaos control to the irreg- 
ular motion of an electron under the combined influence of a 
Coulomb and a magnetic field, the so-called "diamagnetic Ke- 
pler problem" (DKP). We show how to stabilize the classical 
chaotic orbits into regular motion and present a new method 
to follow the unstable orbits as the energy is increased from 
the regular to the chaotic regime. 

PACS numbers: 05.45-a, 05.45. Gg 



The dynamics of an hydrogen atom in a strong mag- 
netic field is known as the "diamagnetic Kepler prob- 
lem" (DKP). After more than two decades of research, 
the DKP is currently the most intensively studied au- 
tonomous, classically chaotic, atomic physics system. As 
first realized by Edmonds the DKP is important in 
such diverse areas as solid state physics, astrophysics, 
and Rydberg atoms and occupies central stage in classi- 
cal and quantum chaos research The classical flow 
covers a wide range of Hamiltonian dynamics reaching 
from bound, nearly integrable behavior to completely 
chaotic and unbound motion as the scaled energy is var- 
ied 1^ . This letter reports the numerical implementation 
of a control scheme to stabilize this chaotic behavior. 
The methodology is sufficiently general to apply to all 
Hamiltonian flows and , with slight modifications, is also 
applicable to ballistic dynamics (e.g. billiards) and area- 
preserving mappings. 

Ever since Ott, Grebogi and York (OGY) Q| intro- 
duced the notion of chaos control a decade ago, the num- 
ber of experimental verifications has continuously and 
rapidly expanded to include a wide variety of applica- 
tions ranging from magnetoelastic ribbons |p|, electronic 
circuits Q , lasers 
convection loops 
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, chemical reactions H to thermal 
The OGY algorithm was written 
explicitly for dissipative systems and most studies (ex- 
perimental and theoretical) have focused their attention 
on this class of systems. Non-trivial modifications are 
needed to adapt the OGY technique to conservative flows 
(or maps) and Lai, Ding and Grebogi [|lO| have provided 
the necessary changes in the context of area-preserving 
mappings. 

Our system is a continuous, 2 degrees of freedom 
(4D phase space) Hamiltonian flow representing the mo- 
tion of an electron under the combined influence of a 
Coulomb and an external magnetic field. We use scaled 
semi-parabolic coordinates and write the resulting scaled 
(pseudo-) Hamiltonian as (for angular momentum L — 0) 



(1) 



The scaled energy e is related to the physical energy E 
by e = 7q ^^"^ E where the parameter 70 = B / Bq denotes 
the strength of the magnetic field relative to the unit 
Bq ~ 2.35 10^ T. The structure of the dynamics of the 
Hamiltonian (Q) depends solely on the value of the scaled 
energy e. For instance, it is known that for —0.5 < e < 
—0.13 the system exhibits bounded motion with mixed 
chaotic and regular motion and for —0.13 < e < 0.0 the 
last large stable island disappears and the dynamics is 
mostly chaotic p^ . 

The dimension reduction (from 4D to 2D) and dis- 
cretization is performed by observing the dynamics on 
the Poincare section defined by /i = 0, /i > 0. The energy 
shell is then mapped to an area bounded by the condi- 
tion pf, — 2e i/^ — A which represents an ellipse in the 
{v,Pv) plane. Figure (1) shows the collection of points 
{^mPv,n} obtained by numerical integration of the equa- 
tions of motion for e = —0.5. One notices, for this energy, 
that the mixed phase space is still mostly regular with a 
small stochastic area filled by the successive piercings of 
one chaotic trajectory. We recall that a periodic trajec- 
tory would show up on the Poincare section as a finite 
number of crossings. 

In attempting to bring order to the DKP dynamics, we 
had to overcome a number of difficulties not encountered 
in previous studies. We had in mind the adaptive control 
of chaotic DKP trajectories under changing conditions 
(e.g. drift in the magnetic field) through the stabilization 
of preselected unstable periodic orbits (UPOs) (moni- 
tored on the Poincare section) via small programmed per- 
turbations to the scaled energy. The difficulties are five- 
fold. First, a typical trajectory spends a lot of time away 
from the Poincare section where the flow is discretized 
and because of the sensitivity of the chaotic dynamics 
we had to device an efficient variable step symplectic in- 
tegrator thereby preserving the geometrical structure of 
the Hamiltonian |l^. Second, the locations of a number 
of UPOs (our target states) and the corresponding Jaco- 
bian matrices for each member of the periodic orbits must 
be obtained numerically. It was found necessary to inter- 
vene at every crossing of the Poincare section to assure a 
robust stabilization. Third, because of anomalous trans- 
port in Hamiltonian phase space p4|, the waiting time 
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for a trajectory to enter the neighborhood of a predefined 
target UPO is unduly long and a targeting strategy is in 
general necessary to steer the trajectory towards the ap- 
propriate region of phase space. This problem has been 
addressed in and is assumed to be solved. Fourth, the 
eigenvalues of area-preserving Jacobians are often com- 
plex and the stable and unstable manifolds, at the core of 
the OGY procedure, are no longer along the directions 
of their eigenvectors. A new method had to be imple- 
mented along the lines described in ]lO| . Fifth, a tracking 
algorithm for updating the control parameters (positions 
of UPOs, Jacobian matrices ...) was designed to insure 
that control is kept while external conditions are chang- 
ing. Part of our solutions to these problems are discussed 
below and the details will be reported elsewhere. 

To locate the UPOs on the Poincare section, we 
have used two methods: a modified recurrence method 
(MRM) and an extension to conservative mappings of a 
stability transform algorithm (STA) introduced in 
The MRM uses the basic idea proposed by Auerbach et 
al. [0 for extracting periodic orbits from a chaotic time 
series {Xiji^i^^f. However, instead of scanning a very 
long ^ 1 time series for pairs of points separated in 
time by m iterates (for a period-m UPO) and within a 
preassigned spatial distance r from each other, M points 
are launched on the Poincare section and integrated for- 
ward until the m-th return on the section where the pair 
of points are compared. In both methods, the accepted 
pairs, a distance r apart, are distributed among the pos- 
sibly different UPO candidates, the center of mass of each 
subset is calculated and the coordinates are stored. For 
a given relative accuracy, 10~'' — 10~^, we find M <^ N 
by several orders of magnitude. In Hamiltonian cases, A^ 
must be at least 10^ in order to reach an acceptable pre- 
cision. In the MRM, one can also easily take advantage 
of the symmetries on the Poincare section for the flow. 
For example, figure (^) shows that the mapping is sym- 
metric with respect to both v and Pu axes. Therefore one 
only needs to spread the initial conditions over one quar- 
ter of phase space for a smaller set of initial conditions 
or a denser distribution. This optimization is essential if 
one considers that between every two consecutive points 
on the Poincare section, a segment of trajectory in phase 
space must be integrated numerically. Furthermore, for 
the DKP, this segment grows rapidly as the scaled energy 
approaches zero where the trajectory spends an increas- 
ingly long time in the branches (they grow as |e|~^^^) of 
the pseudo-potential. 

The STA amounts to the following. Consider a discrete 
chaotic dynamical system S defined in a d-dimensional 
space by 

X„+i=F(X„,p) (2) 

where p is some external parameter. The goal is to con- 
struct from equation (^ other dynamical systems Sk with 



the same periodic orbits (in number and positions), but 
whose stability has changed, i.e unstable orbits become 
stable, and stable orbits remain stable. The new systems 
are defined by the linear transformations (written for a 
period m point, F(™)(X*,p) = X*) 

X„+i = X„ + A Cfc [f(")(X„,p) - X„] (3) 

where < A <C 1 is adjusted to improve convergence and 
Cfc are non-singular constant dx d matrices (see for 
their explicit forms) . Because of the stability of the peri- 
odic orbits of the constructed system Sk, every trajectory 
of Sk converges after a finite number of iterations to a 
periodic m point X*. Per construction, the X* are also 
periodic points of the original system S. We have tested 
this procedure for a number of area preserving mappings 
and flows and the global convergence of the method is ex- 
cellent. For flows, where the mapping in (^) is replaced 
by a numerical integration for m returns to the Poincare 
section, the STA has proved equally reliable. In this re- 
spect, our adaptive symplectic integrator has insured 
that the (long) trajectories away from the Poincare sec- 
tion were kept sufficiently accurate. The position of the 
fixed point in figure (2) (bottom panel) has been reached 
after only 250 iterations starting from an arbitrary initial 
point: this is representative of the efficiency of the STA. 
For the same relative accuracy (10^^), the MRM needs 
approximately 10^ initial points. 

In order to stabilize a chaotic trajectory around one of 
these UPOs, we have implemented a numerical version of 
the OGY method as modified in jl^] for area-preserving 
mappings. This method is believed to be dynamically 
optimal in that it explicitly uses the local geometry of 
the underlying system. Given a dynamical system of the 
type (^), and a target UPO of period to, {X(z,po)}i=i,m 
, at some nominal parameter value poj one characterizes 
the local stable and unstable manifolds by the vectors 
Gg^i and Bu^i respectively as well as their contravariant 
counterparts fg.i and f„.i satisfying iu,i'^u,i — ^s.i'^s.i = 1 
a,nAiu,i'^s,i = ^s,i'^u,i = 0. The stabilizing perturbations 
5pn = Pn — Po are then obtained by firstly linearizing 
the dynamics in a (5-neighborhood of a member of the 
periodic orbit, say X(/c,po) i and around po, namely 

X„+i-X(fc + l,p„)^Ufc [X„-X(fc,p„)] (4) 

where the d x d Jacobian matrix U = Z3x F(X,p) is 
evaluated at [X = X(A:,po)j P = Po] and it is understood 
that ||X„ - X(fc,po)|| < (5 < 1 and \6pn\ < bo|- 

Secondly, the control criterion is imposed that X„+i = 
F(X„,p„) should lie along the stable direction at X(fc + 
l,Po), i-e. ft.,fc+i • [X„+i - X(fc + 1,pq)] = which, 
together with the parametric variation of the periodic 
points, X(/c,po + Sp) X(fc,po) + gfc ^P j leads to the 
following expression for the parameter perturbation at 
the 71-th iteration: 
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The individual Jacobian matrices U,; were obtained 
by least-square minimization ||l8|| of the linear rela- 
tionship between subsequent pairs of iterates on the 
Poincare section for over 1000 points drawn from the 
close neighborhood of each X(i,po)- For the dissipative 
cases, the eigenvectors of the Jacobian matrix U^,™'' = 

DxF(™)(Xfe,po) = I\"=o^ '^k+j would determine the 
manifold directions, but as first realized in pO{ , the con- 
servative cases lead generically to complex eigenvalues 
and eigenvectors and a different strategy is required. The 
solution adopted to generate the unstable (stable) direc- 
tion e„ (Gs) is to rotate an arbitrary unit vector Uq (sq) 
according to u[.'^™-'uo and U^~''™-*So which converge re- 
spectively to Bu^k and Bs^k for g ^ 1. The vectors are 
periodically renormalized to prevent overflow and in gen- 
eral, g < 10 is sufficient for convergence. This procedure 
is very reliable and robust. 

Putting all the ingredients together, we have achieved 
the stabilization of the chaotic behavior of the 
diamagnetic-Kepler Hamiltonian (|^) for different values 
of e and different UPOs. Figure (2) presents a sample of 
results for the control of period-2, period-3 and period- 
1 at e = -0.3, -0.2 and -0.1 respectively. The control 
mechanism was switched on when an iterate entered a 
neighborhood of size 5 = 10~^ ( in the 3 panels, one has 
removed the initial transient) and the perturbations on 
the scaled energy (alternatively the magnetic field) were 
never larger than Se = 10~^. After 3000 iterates, we re- 
lease the control, (5e„ = , and the trajectory returns to 
its natural chaotic behavior. One also notes (left panels) 
that the dynamics is increasingly chaotic as e approaches 
zero and by e = —0.1 no obvious trace of regularity is 
present. We emphasize that in performing the present 
control task no a priori knowledge of the governing map 
on the Poincare section is known and the complete pro- 
cedure is purely numerical from beginning to end. 

Under parameter variation, the location of the UPOs 
changes as well as their associated Jacobian matrices, 
and unless action is taken the control over the chaotic 
dynamics may be lost. The strategy of detecting param- 
eter drift and updating accordingly the control param- 
eters is known as tracking ||l9[| . We have developed a 
new tracking algorithm which is simple and efficient. For 
concreteness, we will concentrate the presentation on the 
updating localization of the UPOs. 

The method is based on a judicious construction of a 
family of extrapolating rational functions 1 20 1 . One starts 
from a number s of UPO coordinates (the seeds), e.g. 
{i^i, €i}i=i^s: from which a rational function, i?o(s,e) is 
constructed. Rq is then evaluate at Cg+i = -I- Ae with 
Ae <C 1 to give an extrapolated value Vs+i- A new set 
of seeds is formed {z^i+i, e^+iji^i^s , removing {z^i,ei} 
and keeping the number of seeds fixed , and a new ex- 



trapolating function i?i(s,e) is calculated. The process 
is repeated until a final value e/ is reached. The same 
procedure may also be applied to the pairs {pu,i, £i}i=i.,s 
to provide initial conditions {vi^p^^i} for the integration 
of the complete phase-space trajectory. Two remarks are 
in order. First, this approach is much more accurate at 
6/ than the direct extrapolation i?o(s,e/), although all 
the original seeds are replaced by extrapolated values af- 
ter s iterations. Second, the algorithm permits to cover 
a wide range of parameter variations beyond the initial 
seeds (3-5 in practice) while keeping high accuracy. Since 
the individual extrapolation steps are chosen small, we 
call the method adiabatic localization (AL) or more gen- 
erally adiabatic tracking (AT). 

We have applied AT to the dissipative Henon map for 
the same range chosen in 119 1. The results are in complete 
agreement with those of and with the precise values 
of the position of UPOs obtained by STA. We have also 
applied this algorithm to the standard map, the Henon- 
Heiles Hamiltonian and the DKP and more complete re- 
sults will be reported at a later time. Figure (^) shows 
the tracking of a period- 1 orbit for the DKP. We started 
from five seeds for e: - 0.3, - 0.29, - 0.28, -0.27, -0.26 
and set Ae — 0.001. The tracking process is extended to 
e — —0.075. The results from the STA confirms the ac- 
curacy of the adiabatic tracking. The right panel clearly 
indicates that the result from AT for e = —0.075 is indeed 
a period 1 orbit. 

We have presented for the first time the application 
of a numerical implementation of the OGY control to a 
realistic conservative flow, the DKP Hamiltonian. On 
our way to a complete solution, we have i. developed 
a new variable step symplectic integrator, ii. modified 
the standard recurrence method and extended a powerful 
technique (STA) to locate the positions of the unstable 
periodic orbits, iii. obtained accurate numerical Jaco- 
bian matrices and the local manifolds, and iv. provided 
an efhcient approach (AT) for tracking the UPOs under 
parameter changes. We are presently considering two im- 
portant mesoscopic applications: nonlinear dynamics in 
electromagnetic tra ps and chaotic ray dynamics in 
lasing microcavities ||22|] . It still remains an open question 
however if manipulations of an external parameter (e.g. 
a magnetic field) to induce stabilization of a classical un- 
stable orbit can be extended to the semi-classical regime, 
for example, for the control of Rydberg wave packet dy- 
namics. Research along these lines is actively being pur- 
sued. 
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FIG. 1. Poincare section of the DKP for e = -0.5. 
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FIG. 2. OGY Control of Diamagnetic Kepler orbits for 

3 scaled energies eo = —0.3, —0.2, —0.1 and periods 2,3, 1 in 
top, middle ,and bottom panncls: {left) Poincare section (PS) 
= 0, /j > showing one chaotic trajectory (filled space) and 
the controlled UPO (black dots); (middle) the stabilized pi, 
or ;/ variable for the first 3000 intersections with the PS be- 
fore control is turned off; (right) corresponding 3D stabilized 
trajectory. 




FIG. 3. (middle) Tracking of an unstable fixed point for - 
0.3< e <- 0.075. Solid line: AT. •: STA. (left) 3D trajectory 
for € = -0.3. (right) 3D trajectory for e = -0.075. The ini- 
tial condition for this trajectory is that obtained by adiabatic 
tracking, i.e. (ly* = -1.46942 , pt = -0.000133596). 
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